#################################################################################
### Racialized Misinformation, Factual Corrections, and Prejudicial Attitudes ###
### [Main Replication Code for Study 2]                                       ###
### Authors: Eddy S. F. Yeung, Joseph Glasgow                                 ###
### Date: June 26, 2025                                                       ###
#################################################################################

### Set-up ----
## Clean the working environment and set the working directory
rm(list = ls())
setwd("~/Desktop/racialized_misinfo/replication/Study 2") # set your working directory here, which should also contain the survey data ("study2_dataset.csv")

## Import the dataset
df <- read.csv("study2_dataset.csv")

## Load the required packages
library(tidyverse)
library(texreg)
library(estimatr)
library(cowplot)
library(grid)
library(gridExtra)
library(interflex)
library(multcomp)
library(broom)
library(ivreg)
library(modelsummary)
library(effectsize)

### Code demographic variables ----
## Age
df$age <- df$yob + 2023 - 2006

## Female (= 1)
df$female <- ifelse(df$sex == 2, 1, 0)

## White (= 1)
df$white <- ifelse(as.integer(grepl("1", df$race)) == 1, 1, 0)

## Married (= 1)
df$married <- ifelse(df$marital == 2, 1, 0)

## Education (less than high school = 0; advanced degree = 5)
df$education <- df$educ - 1

## Income (less than $10,000 = 0; $150,000 or more = 12)
df$income <- df$income - 1
df$income <- ifelse(df$income == 13, 12, df$income)

## Party identification (0 = strong Democrat; 6 = strong Republican)
df <- df %>% mutate(pid = case_when(
  pid_1 == 2 & pid_2d == 1 ~ 0,
  pid_1 == 2 & pid_2d == 2 ~ 1,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 2 ~ 2,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 3 ~ 3,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 1 ~ 4,
  pid_1 == 1 & pid_2r == 2 ~ 5,
  pid_1 == 1 & pid_2r == 1 ~ 6
))

## Democrat (= 1)
df$dem <- ifelse(df$pid == 0 | df$pid == 1 | df$pid == 2, 1, 0)

## Republican (= 1)
df$gop <- ifelse(df$pid == 4 | df$pid == 5 | df$pid == 6, 1, 0)

## Self-reported ideology (0 = extremely liberal; 6 = extremely conservative)
df$ideo <- df$ideo - 1

## Attentiveness (0 = lowest; 2 = highest)
df$screener_2_correct <- ifelse(df$screener_2 == 3, 1, 0)
df$screener_4_correct <- ifelse(df$screener_4 == 2, 1, 0)
df$attentive <- df$screener_2_correct + df$screener_4_correct
table(df$attentive)

### Code outcome variables ----
## Racial resentment (mean effects index with imputation)
rr_matrix <- cbind(df$racial_resent_1, df$racial_resent_2, df$racial_resent_3, 
                   df$racial_resent_4, df$racial_resent_5)
cal_MEI <- 
  function(Z, outcome_mat, to_reorient, reorient = F, greedy = T, impute = F) {
    if(impute == T) {
      R <- 1 * is.na(outcome_mat)
      means_for_imputation <- 
        rbind(apply(outcome_mat[Z == 0, ], MAR = 2, FUN = mean, na.rm = T),
              apply(outcome_mat[Z == 1, ], MAR = 2, FUN = mean, na.rm = T))
      to_impute <- R * means_for_imputation[Z + 1, ]
      outcome_mat[is.na(outcome_mat)] <- 0
      outcome_mat <- outcome_mat + to_impute
    }
    c_mean <- apply(X = outcome_mat[Z == 0, ], MARGIN = 2, FUN = mean, na.rm = T)
    c_sd <- apply(X = outcome_mat[Z == 0, ], MARGIN = 2, FUN = sd, na.rm = T)
    z_score <- t(t(sweep(outcome_mat, 2, c_mean)) / c_sd)
    index_numerator <- rowSums(z_score)
    if(greedy == T) {
      n_outcomes <- rowSums(!is.na(z_score))
    }
    else if(greedy == F){
      n_outcomes <- ncol(outcome_mat)
    }
    index <- index_numerator / n_outcomes
    index <- (index - mean(index[Z == 0], na.rm = T)) / sd(index[Z == 0], na.rm = T)
    return(index)
  }
df <- df %>% 
  mutate(resent = cal_MEI(Z = df$treatment, outcome_mat = rr_matrix, impute = T))

## Explicit prejudice (difference between Blacks' laziness and Whites' laziness)
df$explicit_prej <- df$explicit_prej_2 - df$explicit_prej_1
table(df$explicit_prej)

## Support for welfare (0 = strongly against; 6 = strongly in favor)
df <- df %>% mutate(welfare_support = case_when(
  welfare_support_2b == 1 ~ 0,
  welfare_support_2b == 2 ~ 1,
  welfare_support_2c == 2 ~ 2,
  welfare_support_2c == 3 ~ 3,
  welfare_support_2c == 1 ~ 4,
  welfare_support_2a == 2 ~ 5,
  welfare_support_2a == 1 ~ 6
))
table(df$welfare_support)

## Posttreatment belief of the share of welfare recipients (1 = correct)
df$mp_check <- ifelse(df$fmc == 2, 1, 0)

## Posttreatment belief about the percentage of Whites in the US population
df$pop_white <- df$population_1

## Posttreatment belief about the percentage of Blacks in the US population
df$pop_black <- df$population_2

## Posttreatment belief about the proportion of welfare recipients among Blacks
df$welfare_black <- df$welfare_black_1

## Overestimation of the share of Black recipients of welfare (truth = 21%)
df$overest <- df$prior_2 - 21

df$treatment <- as.factor(df$treatment)

### Figure 3 ----
## Average level of misperception in the full sample
mean(df$prior_2, na.rm = T)
median(df$prior_2, na.rm = T)
length(na.omit(df$prior_2))

## Pretreatment estimate of the proportion of Black welfare recipients by party ID
group_mean_temp1 <- df %>% 
  filter(gop == 0 | gop == 1) %>% 
  group_by(gop) %>% 
  summarize(xvalue = mean(prior_2, na.rm = T))
p0.1 <- 
  ggplot(data = subset(df, gop == 0 | gop == 1), 
         aes(x = prior_2, group = as.factor(gop), fill = as.factor(gop))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp1,
             aes(xintercept = xvalue, color = as.factor(gop)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 21, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated % of Black Recipients") +
  ylab("Density") +
  ggtitle("Prior Beliefs by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))

## Pretreatment estimate of the proportion of Black welfare recipients by racial identity
group_mean_temp2 <- df %>% 
  filter(white == 0 | white == 1) %>% 
  group_by(white) %>% 
  summarize(xvalue = mean(prior_2, na.rm = T))
p0.2 <- 
  ggplot(data = subset(df, white == 0 | white == 1), 
         aes(x = prior_2, group = as.factor(white), fill = as.factor(white))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp2,
             aes(xintercept = xvalue, color = as.factor(white)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 21, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated % of Black Recipients") +
  ylab("Density") +
  ggtitle("Prior Beliefs by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
prior <- plot_grid(p0.1, p0.2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure 3.pdf", prior, width = 8, height = 3)
# ggsave(file = "Figure 3.tif", prior, dpi = 300, width = 8, height = 3)

### Figure 4 and Figure S13 ----
## Explicit prejudice
# Estimate the ATE
prej_ATE <- lm_robust(explicit_prej ~ treatment + age + female + white + married +
                        education + income + ideo + pid,
                      data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(explicit_prej ~ treatment, data = df)

# Compute the general distributions
prej_by_group <- df %>% 
  filter(explicit_prej >= -6 & explicit_prej <= 6) %>% 
  group_by(treatment, explicit_prej) %>% 
  summarize(count = n()) %>% 
  mutate(freq = formattable::percent(count / sum(count)))

# Visualize the results
group_mean_temp1 <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(explicit_prej, na.rm = T))
p1 <- 
  ggplot(data = prej_by_group, 
         aes(x = explicit_prej, y = freq * 100, group = treatment, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(.9), 
           color = "black", alpha = 0.25) +
  geom_vline(data = group_mean_temp1,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Explicit Prejudice against Blacks\n(-6 = lowest; 6 = highest)") +
  ylab("Percentage of Respondents (%)") +
  ggtitle("Explicit Prejudice in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p1 <- p1 + 
  annotate(geom = "text", x = -3, y = 37, 
              label = paste("beta == ", round(prej_ATE$coefficients[2], 2)),
              size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = -3, y = 34.4, 
           label = paste("(SE == ", round(prej_ATE$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)

# HTE by seven-point party ID
out <- interflex(Y = "explicit_prej", D = "treatment", X = "pid",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo"), 
                 data = df, na.rm = TRUE,
                 estimator = "binning", vcov.type = "robust")
p1_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Party Identification\n(0 = strong Democrat; 6 = strong Republican)", 
       ylab = "Marginal Effect on Explicit Prejudice",
       ylim = c(-1.5, 1.5)) +
  ggtitle("HTE on Explicit Prejudice by Party Identification") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))

## Racial resentment
# Estimate the ATE
resent_ATE <- lm_robust(resent ~ treatment + age + female + white + married +
                          education + income + ideo + pid,
                        data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(resent ~ treatment, data = df)

# Visualize the results
group_mean_temp2 <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(resent, na.rm = T))
p2 <- 
  ggplot(data = df, aes(x = resent, group = treatment, fill = treatment)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp2,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Racial Resentment\n(Mean Effects Index)") +
  ylab("Density") +
  xlim(-2.1, 2.1) +
  ggtitle("Racial Resentment in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p2 <- p2 + 
  annotate(geom = "text", x = -1, y = 0.40, 
           label = paste("beta == ", round(resent_ATE$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = -1, y = 0.372, 
           label = paste("(SE == ", round(resent_ATE$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
main_ATE <- plot_grid(p1, p2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure 4.pdf", main_ATE, width = 9, height = 4.5)
# ggsave(file = "Figure 4.tif", main_ATE, dpi = 300, width = 9, height = 4.5)

# HTE by seven-point party ID
out <- interflex(Y = "resent", D = "treatment", X = "pid",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo"), 
                 data = df, na.rm = TRUE, 
                 estimator = "binning", vcov.type = "robust")
p2_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Party Identification\n(0 = strong Democrat; 6 = strong Republican)", 
       ylab = "Marginal Effect on Racial Resentment",
       ylim = c(-.75, .75)) +
  ggtitle("HTE on Racial Resentment by Party Identification") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))
HTE_pid <- plot_grid(p1_HTE, p2_HTE, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S13.pdf", HTE_pid, width = 9, height = 4.5)

### Figure S23 ----
## Support for welfare
# Estimate the ATE
ATE_welfare_support <- lm_robust(welfare_support ~ treatment + age + female + white + 
                                   married + education + income + ideo + pid,
                                 data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(welfare_support ~ treatment, data = df)

# Visualize the ATE
summary_welfare_support <- df %>% 
  group_by(treatment) %>% 
  do(tidy(lm_robust(welfare_support ~ 1, data = .))) %>% 
  mutate(welfare_support = estimate)
p <- 
  ggplot(summary_welfare_support, aes(x = treatment, y = welfare_support)) +
  geom_point(data = df, size = 1, alpha = .15, color = "grey",
             position = position_jitter(width = .1, height = .2, seed = 1)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), linewidth = .5, width = 0) +
  theme_classic() +
  scale_x_discrete(labels = c("0" = "Control", "1" = "Treatment")) +
  xlab("") +
  ylab("Support for Welfare\n(0 = lowest; 6 = highest)") +
  coord_cartesian(ylim = c(-.25, 6.25)) +
  theme(text = element_text(family = "Times", size = 12),
        axis.text = element_text(family = "Times", size = 12))
p <- p + 
  annotate(geom = "text", x = 1.5, y = 5.15, 
           label = paste("beta == ", round(ATE_welfare_support$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 1.5, y = 4.85, 
           label = paste("(SE == ", round(ATE_welfare_support$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
ggsave(file = "Figure S23.pdf", width = 5, height = 5)

### Table S6 ----
## Identify respondents who believed that Black welfare recipients outnumbered White welfare recipients
df$misper <- ifelse(df$prior_2 > df$prior_1, 1, 0)

## Analyze the correlations between misperceptions and explicit prejudice / racial resentment
mod1.1 <- lm_robust(explicit_prej ~ treatment, 
                    data = df)
mod1.2 <- lm_robust(explicit_prej ~ treatment + age + female + white + 
                      married + education + income + ideo + pid, 
                    data = df)
mod1.3 <- lm_robust(explicit_prej ~ treatment + age + female + white + 
                      married + education + income + ideo + pid + misper, 
                    data = df)
mod1.4 <- lm_robust(explicit_prej ~ treatment + age + female + white + 
                      married + education + income + ideo + pid + misper, 
                    data = df, subset = attentive == 2)
mod2.1 <- lm_robust(resent ~ treatment, 
                    data = df)
mod2.2 <- lm_robust(resent ~ treatment + age + female + white + 
                      married +education + income + ideo + pid, 
                    data = df)
mod2.3 <- lm_robust(resent ~ treatment + age + female + white + 
                      married +education + income + ideo + pid + misper, 
                    data = df)
mod2.4 <- lm_robust(resent ~ treatment + age + female + white + 
                      married +education + income + ideo + pid + misper, 
                    data = df, subset = attentive == 2)

## Tabularize the results
texreg(list(mod1.1, mod1.2, mod1.3, mod1.4, mod2.1, mod2.2, mod2.3, mod2.4),
       custom.coef.names = c("Intercept", "Treatment (= 1)", "Age", 
                             "Female (= 1)", "White (= 1)", "Married (= 1)", 
                             "Education (5 = highest)", "Income (12 = highest)", 
                             "Self-Reported Ideology", "Party Identification",
                             "Misperception (= 1)"),
       custom.header = list("Dependent Variable" = 1:8),
       custom.model.names = c("Prejudice", "Prejudice", "Prejudice", "Prejudice",
                              "Resentment", "Resentment", "Resentment", "Resentment"),
       stars = numeric(0),
       include.adjrs = TRUE, include.ci = FALSE, include.rmse = FALSE, 
       caption = "Misperceptions about Welfare Recipiency, Factual Corrections, and Racial Prejudice")

### Figure S14 ----
## Standardize the outcome variables
df <- df %>% mutate(
  explicit_prej_z = scale(explicit_prej),
  resent_z = scale(resent),
  welfare_support_z = scale(welfare_support)
)

## Create an empty data frame to store the estimates
df_race <- as.data.frame(matrix(NA, nrow = 6, ncol = 5))
df_race <- df_race %>% 
  rename(race = V1, dv = V2, cate = V3, lower_ci = V4, upper_ci = V5)
df_race$race <- c("Nonwhite", "White", "Nonwhite", "White", "Nonwhite", "White")
df_race$dv <- c(1, 1, 2, 2, 3, 3)

## DV = explicit prejudice
temp1 <- lm_robust(explicit_prej_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp1_nonwhite <- glht(temp1, linfct = c("treatment1 = 0"))
CATE_temp1_white <- glht(temp1, linfct = c("treatment1 + treatment1:white = 0"))
df_race[1, 3:5] <- confint(CATE_temp1_nonwhite)$confint
df_race[2, 3:5] <- confint(CATE_temp1_white)$confint

## DV = racial resentment
temp2 <- lm_robust(resent_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp2_nonwhite <- glht(temp2, linfct = c("treatment1 = 0"))
CATE_temp2_white <- glht(temp2, linfct = c("treatment1 + treatment1:white = 0"))
df_race[3, 3:5] <- confint(CATE_temp2_nonwhite)$confint
df_race[4, 3:5] <- confint(CATE_temp2_white)$confint

## DV = welfare support
temp3 <- lm_robust(welfare_support_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp3_nonwhite <- glht(temp3, linfct = c("treatment1 = 0"))
CATE_temp3_white <- glht(temp3, linfct = c("treatment1 + treatment1:white = 0"))
df_race[5, 3:5] <- confint(CATE_temp3_nonwhite)$confint
df_race[6, 3:5] <- confint(CATE_temp3_white)$confint

## Visualize the CATEs by racial identity and by dependent variable (all results)
df_race$dv <- factor(df_race$dv,
                     levels = c(1:3),
                     labels = c("Explicit Prejudice", "Racial Resentment", "Welfare Support"))
p <- 
  ggplot(data = df_race, aes(x = dv, y = cate, color = race, shape = race)) +
  geom_point(position = position_dodge(.25), size = 3) +
  scale_color_manual(values = c("grey25", "grey75")) +
  scale_shape_manual(values = c(16, 17)) +
  geom_errorbar(width = 0, aes(ymin = lower_ci, ymax = upper_ci), 
                position = position_dodge(.25)) +
  xlab("") + 
  ylab("Conditional Average Treatment Effect") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.99, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-.5, .5))
ggsave(file = "Figure S14.pdf", width = 8, height = 4)

### Figure S6 ----
## % of respondents disagreeing that "Black recipients > White recipients" in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(mp_check, na.rm = T))
lm_robust(mp_check ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## % of respondents disagreeing that "Black recipients > White recipients" by party ID
df_manip <- as.data.frame(matrix(NA, nrow = 4, ncol = 5))
df_manip <- df_manip %>% 
  rename(gop = V1, group = V2, percent = V3, lower_ci = V4, upper_ci = V5)
df_manip$gop <- c("Republican", "Republican", "Non-Republican", "Non-Republican")
df_manip$group <- c("Control", "Treatment", "Control", "Treatment")
temp1 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$gop == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$gop == 1])),
                   correct = F)
temp2 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$gop == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$gop == 1])),
                   correct = F)
temp3 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$gop == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$gop == 0])),
                   correct = F)
temp4 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$gop == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$gop == 0])),
                   correct = F)
df_manip[1, 3] <- temp1$estimate
df_manip[2, 3] <- temp2$estimate
df_manip[3, 3] <- temp3$estimate
df_manip[4, 3] <- temp4$estimate
df_manip[1, 4:5] <- temp1$conf.int
df_manip[2, 4:5] <- temp2$conf.int
df_manip[3, 4:5] <- temp3$conf.int
df_manip[4, 4:5] <- temp4$conf.int
df_manip$gop <- factor(df_manip$gop, levels = c("Republican", "Non-Republican"))
p_manip_pid <- 
  ggplot(data = df_manip, aes(x = gop, y = percent * 100, fill = group)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = lower_ci * 100, ymax = upper_ci * 100), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Percent of Respondents Disagreeing that\nBlack Welfare Recipients > White Welfare Recipients") +
  ggtitle("Direct Manipulation Check by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))

## % of respondents disagreeing that "Black recipients > White recipients" by racial identity
df_manip <- as.data.frame(matrix(NA, nrow = 4, ncol = 5))
df_manip <- df_manip %>% 
  rename(gop = V1, group = V2, percent = V3, lower_ci = V4, upper_ci = V5)
df_manip$gop <- c("White", "White", "Nonwhite", "Nonwhite")
df_manip$group <- c("Control", "Treatment", "Control", "Treatment")
temp1 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$white == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$white == 1])),
                   correct = F)
temp2 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$white == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$white == 1])),
                   correct = F)
temp3 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$white == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$white == 0])),
                   correct = F)
temp4 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$white == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$white == 0])),
                   correct = F)
df_manip[1, 3] <- temp1$estimate
df_manip[2, 3] <- temp2$estimate
df_manip[3, 3] <- temp3$estimate
df_manip[4, 3] <- temp4$estimate
df_manip[1, 4:5] <- temp1$conf.int
df_manip[2, 4:5] <- temp2$conf.int
df_manip[3, 4:5] <- temp3$conf.int
df_manip[4, 4:5] <- temp4$conf.int
df_manip$gop <- factor(df_manip$gop, levels = c("White", "Nonwhite"))
p_manip_race <- 
  ggplot(data = df_manip, aes(x = gop, y = percent * 100, fill = group)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = lower_ci * 100, ymax = upper_ci * 100), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Percent of Respondents Disagreeing that\nBlack Welfare Recipients > White Welfare Recipients") +
  ggtitle("Direct Manipulation Check by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))
manip_check <- plot_grid(p_manip_pid, p_manip_race, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S6.pdf", manip_check, width = 9, height = 4.5)

### Figure S7 ----
## Posttreatment estimate of percent of welfare recipients among Blacks in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(welfare_black, na.rm = T))
lm_robust(welfare_black ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## Posttreatment estimate of percent of welfare recipients among Blacks by party ID
summary_welfare_black <- df %>% 
  group_by(treatment, gop) %>% 
  do(tidy(lm_robust(welfare_black ~ 1, data = .))) %>% 
  mutate(welfare_black = estimate) %>% 
  na.omit()
summary_welfare_black$gop <- factor(summary_welfare_black$gop, levels = c(1, 0), 
                                    labels = c("Republican", "Non-Republican"))
summary_welfare_black$treatment <- factor(summary_welfare_black$treatment, levels = c(0, 1), 
                                          labels = c("Control", "Treatment"))
p_manip_pid_2 <- 
  ggplot(data = summary_welfare_black, aes(x = gop, y = welfare_black, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Estimated % of Welfare Recipients\namong the Black Population") +
  ggtitle("Posterior Beliefs by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))

## Posttreatment estimate of percent of welfare recipients among Blacks by racial identity
summary_welfare_black <- df %>% 
  subset(white == 0 | white == 1) %>% 
  group_by(treatment, white) %>% 
  do(tidy(lm_robust(welfare_black ~ 1, data = .))) %>% 
  mutate(welfare_black = estimate) %>% 
  na.omit()
summary_welfare_black$white <- factor(summary_welfare_black$white, levels = c(1, 0), 
                                      labels = c("White", "Nonwhite"))
summary_welfare_black$treatment <- factor(summary_welfare_black$treatment, levels = c(0, 1), 
                                          labels = c("Control", "Treatment"))
p_manip_race_2 <- 
  ggplot(data = summary_welfare_black, aes(x = white, y = welfare_black, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Estimated % of Welfare Recipients\namong the Black Population") +
  ggtitle("Posterior Beliefs by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))
manip_check <- plot_grid(p_manip_pid_2, p_manip_race_2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S7.pdf", manip_check, width = 9, height = 4.5)

### Figure S24 ----
ATE_pop_black <- lm_robust(pop_black ~ treatment + age + female + white + 
                             married + education + income + ideo + pid,
                           data = df)
group_mean_pop_black <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(pop_black, na.rm = T))
p1 <- 
  ggplot(data = df, aes(x = pop_black, group = treatment, fill = treatment)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_pop_black,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 14.2, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated % of Blacks in the US Population") +
  ylab("Density") +
  ggtitle("Posterior Beliefs about Black Population in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p1 <- p1 + 
  annotate(geom = "text", x = 50, y = 0.032, 
           label = paste("beta == ", round(ATE_pop_black$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 50, y = 0.030, 
           label = paste("(SE == ", round(ATE_pop_black$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
ggsave(file = "Figure S24.pdf", p1, width = 6, height = 4.5)

### Figure S15 ----
## HTE on explicit prejudice
out <- interflex(Y = "explicit_prej", D = "treatment", X = "overest",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                 data = df, na.rm = TRUE,
                 estimator = "kernel", vcov.type = "robust")
p1_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Explicit Prejudice",
       ylim = c(-2, 2)) +
  ggtitle("HTE on Explicit Prejudice by Level of Overestimation") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))

## HTE on racial resentment
out <- interflex(Y = "resent", D = "treatment", X = "overest",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                 data = df, na.rm = TRUE,
                 estimator = "kernel", vcov.type = "robust")
p2_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Racial Resentment",
       ylim = c(-1, 1)) +
  ggtitle("HTE on Racial Resentment by Level of Overestimation") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))
HTE_overest <- plot_grid(p1_HTE, p2_HTE, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S15.pdf", HTE_overest, width = 10, height = 4.5)

### Exploratory analysis: unpacking the ATE on explicit prejudice ----
## ATE on Whites' laziness
lm_robust(explicit_prej_1 ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## ATE on Blacks' laziness
lm_robust(explicit_prej_2 ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

### Table S9 ----
## CACE on explicit prejudice
cace_1 <- ivreg(explicit_prej ~ mp_check | treatment, data = df)
cace_2 <- ivreg(explicit_prej ~ mp_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## CACE on Whites' laziness
cace_3 <- ivreg(explicit_prej_1 ~ mp_check | treatment, data = df)
cace_4 <- ivreg(explicit_prej_1 ~ mp_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## CACE on Blacks' laziness
cace_5 <- ivreg(explicit_prej_2 ~ mp_check | treatment, data = df)
cace_6 <- ivreg(explicit_prej_2 ~ mp_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## CACE on racial resentment
cace_7 <- ivreg(resent ~ mp_check | treatment, data = df)
cace_8 <- ivreg(resent ~ mp_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## Summarize the results
m_list <- list(IV1 = cace_1, IV2 = cace_2, IV3 = cace_3, IV4 = cace_4,
               IV5 = cace_5, IV6 = cace_6, IV7 = cace_7, IV8 = cace_8)
msummary(m_list, output = "markdown")
